execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

ex12

hierarchical model  
class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)
n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
  scale_shape_manual(values=1:k)

qplot(x,y,col=c)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
estimate as no class
y~N(b00+b10*x,s)

ex8-0.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
mdl=cmdstan_model('./ex8-0.stan') 
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)
## Initial log joint probability = -123310 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpFvED78/model-9918b9cc8.stan', line 14, column 2 to column 22) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       28      -161.995   0.000789827    0.00109828           1           1       75    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__    -162.00
##    b0        15.75
##    b1         0.83
##    s         15.49
##    m0[1]     18.81
##    m0[2]     20.27
##    m0[3]     26.60
##    m0[4]     19.26
##    m0[5]     18.67
##    m0[6]     23.39
##    m0[7]     16.82
##    m0[8]     27.36
##    m0[9]     20.53
##    m0[10]    32.04
##    m0[11]    15.97
##    m0[12]    25.89
##    m0[13]    19.34
##    m0[14]    28.25
##    m0[15]    20.16
##    m0[16]    26.65
##    m0[17]    29.80
##    m0[18]    23.55
##    m0[19]    24.53
##    m0[20]    22.42
##    m0[21]    30.11
##    m0[22]    32.20
##    m0[23]    23.57
##    m0[24]    17.91
##    m0[25]    25.88
##    m0[26]    17.44
##    m0[27]    27.70
##    m0[28]    16.13
##    m0[29]    17.18
##    m0[30]    20.10
##    m0[31]    22.97
##    m0[32]    23.13
##    m0[33]    18.60
##    m0[34]    28.39
##    m0[35]    23.61
##    m0[36]    18.16
##    m0[37]    25.54
##    m0[38]    17.69
##    m0[39]    17.97
##    m0[40]    29.07
##    m0[41]    23.31
##    m0[42]    28.32
##    m0[43]    30.32
##    m0[44]    30.30
##    m0[45]    20.38
##    m0[46]    23.12
##    m0[47]    24.49
##    m0[48]    29.73
##    m0[49]    22.21
##    m0[50]    19.29
##    y0[1]     13.63
##    y0[2]     48.27
##    y0[3]      7.22
##    y0[4]     17.49
##    y0[5]     17.68
##    y0[6]     36.86
##    y0[7]      8.77
##    y0[8]     47.78
##    y0[9]     22.35
##    y0[10]    21.51
##    y0[11]     1.30
##    y0[12]    26.93
##    y0[13]    33.22
##    y0[14]    33.51
##    y0[15]    -6.02
##    y0[16]    25.54
##    y0[17]    39.50
##    y0[18]    27.53
##    y0[19]    23.70
##    y0[20]    34.32
##    y0[21]    47.03
##    y0[22]    35.35
##    y0[23]    42.43
##    y0[24]    -4.79
##    y0[25]    -0.08
##    y0[26]    29.59
##    y0[27]    18.88
##    y0[28]     4.11
##    y0[29]    29.39
##    y0[30]    26.19
##    y0[31]     2.79
##    y0[32]    12.46
##    y0[33]    33.39
##    y0[34]    38.35
##    y0[35]    36.89
##    y0[36]     1.64
##    y0[37]    43.46
##    y0[38]     9.08
##    y0[39]    30.30
##    y0[40]    18.68
##    y0[41]    16.52
##    y0[42]     4.32
##    y0[43]    20.40
##    y0[44]    34.81
##    y0[45]    32.87
##    y0[46]     1.62
##    y0[47]    27.84
##    y0[48]    34.31
##    y0[49]    33.94
##    y0[50]     2.92
##    m1[1]     15.97
##    m1[2]     17.77
##    m1[3]     19.58
##    m1[4]     21.38
##    m1[5]     23.18
##    m1[6]     24.98
##    m1[7]     26.79
##    m1[8]     28.59
##    m1[9]     30.39
##    m1[10]    32.20
##    y1[1]      9.86
##    y1[2]      6.20
##    y1[3]     39.43
##    y1[4]      7.12
##    y1[5]     41.92
##    y1[6]      1.27
##    y1[7]     24.37
##    y1[8]     33.12
##    y1[9]     28.90
##    y1[10]     7.75
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -160.73 -160.44  1.19  0.98 -163.01 -159.44 1.00      718      774
##    b0       16.26   16.15  4.31  4.24    9.05   23.14 1.01      312      387
##    b1        0.78    0.79  0.40  0.41    0.12    1.43 1.00      414      684
##    s        16.25   16.08  1.69  1.71   13.83   19.20 1.00     2313     1543
##    m0[1]    19.16   19.14  3.14  3.08   13.83   24.15 1.01      319      433
##    m0[2]    20.53   20.56  2.69  2.68   15.95   24.94 1.01      345      505
##    m0[3]    26.53   26.52  2.73  2.59   22.00   30.99 1.00     1466     1193
##    m0[4]    19.58   19.59  2.99  2.97   14.48   24.35 1.01      324      452
##    m0[5]    19.02   19.00  3.19  3.13   13.62   24.08 1.01      318      433
##    m0[6]    23.49   23.53  2.24  2.21   19.72   27.09 1.00      652     1119
##    m0[7]    17.27   17.21  3.87  3.79   10.75   23.48 1.01      312      387
##    m0[8]    27.25   27.25  2.96  2.82   22.46   32.07 1.00     1393     1209
##    m0[9]    20.78   20.82  2.62  2.61   16.28   25.08 1.01      354      555
##    m0[10]   31.68   31.61  4.75  4.81   24.03   39.48 1.00      794     1233
##    m0[11]   16.46   16.35  4.22  4.15    9.42   23.23 1.01      312      384
##    m0[12]   25.86   25.86  2.55  2.39   21.52   29.96 1.00     1431     1193
##    m0[13]   19.66   19.67  2.96  2.95   14.60   24.41 1.01      325      454
##    m0[14]   28.09   28.08  3.25  3.16   22.83   33.35 1.00     1272     1225
##    m0[15]   20.43   20.46  2.72  2.70   15.78   24.88 1.01      342      497
##    m0[16]   26.57   26.58  2.74  2.60   22.02   31.06 1.00     1463     1193
##    m0[17]   29.56   29.51  3.83  3.79   23.34   35.81 1.00      999     1328
##    m0[18]   23.64   23.69  2.24  2.19   19.91   27.26 1.00      692     1174
##    m0[19]   24.57   24.59  2.31  2.21   20.69   28.34 1.00     1059     1259
##    m0[20]   22.57   22.61  2.28  2.26   18.76   26.30 1.01      488     1023
##    m0[21]   29.86   29.80  3.96  3.93   23.38   36.34 1.00      957     1290
##    m0[22]   31.83   31.76  4.82  4.86   24.10   39.73 1.00      784     1218
##    m0[23]   23.66   23.70  2.24  2.18   19.94   27.29 1.00      697     1174
##    m0[24]   18.30   18.32  3.46  3.39   12.48   23.83 1.01      314      421
##    m0[25]   25.85   25.85  2.55  2.39   21.51   29.95 1.00     1431     1193
##    m0[26]   17.85   17.82  3.63  3.56   11.72   23.68 1.01      313      417
##    m0[27]   27.57   27.55  3.07  2.96   22.61   32.49 1.00     1350     1241
##    m0[28]   16.62   16.51  4.15  4.07    9.70   23.29 1.01      312      381
##    m0[29]   17.61   17.58  3.73  3.66   11.30   23.60 1.01      312      406
##    m0[30]   20.38   20.41  2.74  2.71   15.71   24.85 1.01      341      494
##    m0[31]   23.09   23.15  2.25  2.22   19.31   26.70 1.00      567     1023
##    m0[32]   23.24   23.28  2.24  2.22   19.46   26.84 1.00      597      977
##    m0[33]   18.95   18.94  3.21  3.14   13.51   24.05 1.01      318      437
##    m0[34]   28.22   28.21  3.30  3.20   22.88   33.53 1.00     1249     1225
##    m0[35]   23.70   23.73  2.24  2.18   19.98   27.33 1.00      708     1188
##    m0[36]   18.54   18.54  3.36  3.27   12.88   23.91 1.01      315      425
##    m0[37]   25.52   25.51  2.47  2.32   21.34   29.44 1.00     1438     1281
##    m0[38]   18.09   18.08  3.54  3.48   12.14   23.78 1.01      313      421
##    m0[39]   18.35   18.37  3.44  3.35   12.59   23.86 1.01      314      418
##    m0[40]   28.86   28.83  3.55  3.45   23.16   34.54 1.00     1115     1241
##    m0[41]   23.42   23.46  2.24  2.22   19.64   27.02 1.00      636      998
##    m0[42]   28.16   28.15  3.28  3.18   22.85   33.46 1.00     1261     1225
##    m0[43]   30.05   30.01  4.04  3.99   23.39   36.63 1.00      935     1257
##    m0[44]   30.03   29.99  4.03  3.99   23.39   36.60 1.00      937     1257
##    m0[45]   20.64   20.67  2.66  2.65   16.10   25.02 1.01      349      527
##    m0[46]   23.23   23.27  2.24  2.22   19.45   26.83 1.00      595     1052
##    m0[47]   24.53   24.54  2.30  2.21   20.66   28.29 1.00     1038     1259
##    m0[48]   29.49   29.45  3.81  3.77   23.32   35.69 1.00     1008     1328
##    m0[49]   22.37   22.42  2.31  2.26   18.54   26.19 1.01      463      894
##    m0[50]   19.60   19.62  2.98  2.97   14.50   24.36 1.01      325      454
##    y0[1]    19.25   19.35 17.04 16.54   -9.56   46.90 1.00     1887     1788
##    y0[2]    20.80   21.21 16.42 15.66   -6.19   47.92 1.00     2020     2024
##    y0[3]    26.76   26.86 16.29 16.08    0.12   54.35 1.00     1900     1853
##    y0[4]    19.25   18.96 17.00 16.36   -8.24   48.38 1.00     1719     1808
##    y0[5]    19.08   19.37 16.85 16.54   -9.69   46.30 1.00     1568     1959
##    y0[6]    23.38   23.40 16.11 16.06   -3.52   49.31 1.00     1956     2053
##    y0[7]    17.86   17.91 16.77 16.48  -10.07   45.99 1.00     1669     1796
##    y0[8]    27.67   26.96 16.53 16.14    0.86   56.08 1.00     1791     1682
##    y0[9]    19.74   19.67 16.55 16.00   -8.06   47.63 1.00     1902     1971
##    y0[10]   31.37   31.41 16.98 16.86    3.11   59.53 1.00     1865     1689
##    y0[11]   16.34   15.98 17.13 17.20  -11.82   43.77 1.00     1441     1643
##    y0[12]   26.15   25.84 16.93 16.36   -1.19   54.60 1.00     2027     1918
##    y0[13]   19.24   19.26 16.55 16.60   -7.95   46.74 1.00     1840     1835
##    y0[14]   28.38   28.64 16.36 15.75    2.53   55.89 1.00     1993     1816
##    y0[15]   20.36   20.40 16.78 16.67   -7.53   46.79 1.00     1873     1739
##    y0[16]   26.86   26.52 16.52 16.46   -0.23   54.49 1.00     1892     1892
##    y0[17]   29.82   29.62 16.76 16.15    3.13   57.92 1.00     1954     1821
##    y0[18]   24.07   24.25 16.29 16.11   -3.11   51.05 1.00     1945     1969
##    y0[19]   24.76   24.34 16.60 16.56   -1.91   52.91 1.00     2021     1777
##    y0[20]   22.10   22.33 16.22 16.70   -5.52   48.30 1.00     2014     1953
##    y0[21]   30.35   30.15 16.83 16.59    2.88   58.34 1.00     1904     1955
##    y0[22]   32.44   31.90 17.31 16.27    3.84   60.81 1.00     1691     1836
##    y0[23]   24.16   24.32 16.56 15.90   -3.63   50.73 1.00     1706     1924
##    y0[24]   18.46   18.62 16.68 16.26   -9.47   45.80 1.00     1375     1908
##    y0[25]   25.92   25.64 16.49 16.38   -1.03   53.84 1.00     1863     2011
##    y0[26]   17.71   17.57 16.69 17.02   -9.34   44.45 1.00     1871     1670
##    y0[27]   27.14   26.79 16.54 16.33   -0.34   53.41 1.00     1981     1879
##    y0[28]   16.79   16.71 17.02 17.02  -11.02   44.82 1.00     1684     1907
##    y0[29]   17.85   17.99 16.66 16.43   -9.77   45.43 1.00     1551     1924
##    y0[30]   19.80   20.38 16.46 16.79   -7.80   45.80 1.00     1783     1884
##    y0[31]   23.29   22.86 15.75 15.51   -1.89   48.90 1.00     1948     1743
##    y0[32]   23.40   23.30 16.33 16.39   -3.05   50.38 1.00     2105     1926
##    y0[33]   19.20   19.43 16.38 16.87   -8.45   45.16 1.00     2012     1933
##    y0[34]   27.97   28.18 16.30 15.74    0.62   53.93 1.00     1899     1893
##    y0[35]   23.88   24.19 17.05 16.79   -4.74   50.85 1.00     1778     1837
##    y0[36]   18.12   17.91 17.05 16.85   -9.22   45.99 1.00     1637     1837
##    y0[37]   25.47   25.27 16.72 16.68   -1.86   52.40 1.00     2109     1863
##    y0[38]   18.23   17.82 17.11 17.23   -9.61   45.97 1.00     1572     1563
##    y0[39]   18.54   18.56 16.15 16.07   -7.86   44.96 1.00     1968     1931
##    y0[40]   28.78   28.25 16.73 15.93    2.16   56.99 1.00     1923     1971
##    y0[41]   23.64   23.61 15.97 15.43   -3.08   49.21 1.00     1826     1761
##    y0[42]   28.40   28.59 16.71 17.06    1.11   55.35 1.00     2004     1852
##    y0[43]   29.94   29.64 17.02 16.68    2.75   58.31 1.00     1953     1865
##    y0[44]   29.52   29.71 16.93 16.48    1.69   57.36 1.00     1947     1954
##    y0[45]   20.25   20.39 17.09 16.78   -8.34   48.28 1.00     1870     1947
##    y0[46]   22.50   22.51 16.51 15.75   -4.76   49.41 1.00     2042     1930
##    y0[47]   24.86   25.15 16.93 16.70   -2.41   52.36 1.00     1767     1990
##    y0[48]   29.51   29.10 17.02 16.94    2.16   57.76 1.00     1853     1845
##    y0[49]   22.32   22.30 17.09 17.16   -6.49   50.66 1.00     1795     1980
##    y0[50]   19.13   18.77 16.93 16.45   -7.84   47.99 1.00     1497     1851
##    m1[1]    16.46   16.35  4.22  4.15    9.42   23.23 1.01      312      384
##    m1[2]    18.17   18.17  3.51  3.44   12.27   23.79 1.01      313      421
##    m1[3]    19.88   19.90  2.89  2.87   14.93   24.59 1.01      329      464
##    m1[4]    21.58   21.60  2.43  2.41   17.42   25.66 1.01      395      734
##    m1[5]    23.29   23.33  2.24  2.21   19.53   26.89 1.00      608      999
##    m1[6]    25.00   25.04  2.37  2.26   20.99   28.81 1.00     1271     1321
##    m1[7]    26.71   26.72  2.78  2.63   22.08   31.26 1.00     1451     1193
##    m1[8]    28.41   28.39  3.38  3.28   22.91   33.84 1.00     1208     1258
##    m1[9]    30.12   30.07  4.07  4.04   23.38   36.74 1.00      928     1257
##    m1[10]   31.83   31.76  4.82  4.86   24.10   39.73 1.00      784     1218
##    y1[1]    16.70   17.03 16.87 16.77  -12.66   44.24 1.00     1433     1929
##    y1[2]    18.04   17.95 17.20 16.34  -10.53   46.42 1.00     1989     1919
##    y1[3]    19.96   20.21 16.48 15.70   -7.69   47.03 1.00     1813     1873
##    y1[4]    21.21   20.85 16.10 15.65   -4.89   47.89 1.00     1969     1854
##    y1[5]    23.57   24.16 16.27 16.08   -2.80   49.78 1.00     1687     1587
##    y1[6]    25.29   25.04 16.44 15.65   -0.89   52.51 1.00     1958     1895
##    y1[7]    27.22   27.20 16.34 15.96    1.17   54.29 1.00     1957     1946
##    y1[8]    28.20   28.45 16.68 16.22    2.01   55.57 1.00     2004     1948
##    y1[9]    30.78   30.85 16.60 16.24    3.07   57.94 1.00     1951     2054
##    y1[10]   32.02   31.93 16.69 16.20    5.25   59.00 1.00     1781     1767

estimate as independent class  
all b0l,b1l do not have a distribution  
class c=1~k 
yc~N(b0c+b1c*x,s)  

ex12-1.stan

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  vector[K] b0;
  vector[K] b1;
  real<lower=0> s;
}
model{
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-1.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -881.851 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -94.2744    0.00374625      0.119494      0.8574      0.8574      111    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199      -94.2546    0.00130082    0.00425903           1           1      217    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      200      -94.2546   0.000716441    0.00153958           1           1      218    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -94.25
##    b0[1]     11.44
##    b0[2]      6.49
##    b0[3]      3.28
##    b0[4]     10.54
##    b0[5]      8.73
##    b0[6]     14.00
##    b0[7]     20.87
##    b0[8]      5.91
##    b0[9]     12.52
##    b1[1]      0.41
##    b1[2]      2.27
##    b1[3]      2.18
##    b1[4]      1.96
##    b1[5]      2.09
##    b1[6]      0.23
##    b1[7]      4.00
##    b1[8]      1.50
##    b1[9]     -0.66
##    s          4.00
##    m0[1]     16.48
##    m0[2]     42.74
##    m0[3]     36.22
##    m0[4]     12.53
##    m0[5]     10.17
##    m0[6]     27.51
##    m0[7]     13.07
##    m0[8]      3.19
##    m0[9]     15.31
##    m0[10]    18.48
##    m0[11]     7.09
##    m0[12]    29.99
##    m0[13]    38.25
##    m0[14]     2.47
##    m0[15]     8.97
##    m0[16]    16.99
##    m0[17]    43.80
##    m0[18]    27.95
##    m0[19]    16.41
##    m0[20]    17.99
##    m0[21]    90.42
##    m0[22]    -0.70
##    m0[23]     6.23
##    m0[24]    14.19
##    m0[25]    16.78
##    m0[26]    29.03
##    m0[27]    34.75
##    m0[28]     6.60
##    m0[29]    12.35
##    m0[30]    18.47
##    m0[31]    26.36
##    m0[32]    19.27
##    m0[33]    12.85
##    m0[34]    17.47
##    m0[35]    16.16
##    m0[36]    32.55
##    m0[37]    23.64
##    m0[38]    12.40
##    m0[39]    31.59
##    m0[40]    42.07
##    m0[41]     6.44
##    m0[42]    28.69
##    m0[43]    32.30
##    m0[44]    45.57
##    m0[45]    43.27
##    m0[46]    26.75
##    m0[47]    16.40
##    m0[48]    44.95
##    m0[49]    14.65
##    m0[50]    37.98
##    y0[1]     15.46
##    y0[2]     45.32
##    y0[3]     35.44
##    y0[4]     19.34
##    y0[5]      7.07
##    y0[6]     30.07
##    y0[7]     19.67
##    y0[8]      0.66
##    y0[9]     21.42
##    y0[10]    24.81
##    y0[11]     6.86
##    y0[12]    38.44
##    y0[13]    38.93
##    y0[14]     3.97
##    y0[15]     0.26
##    y0[16]    19.08
##    y0[17]    48.03
##    y0[18]    23.31
##    y0[19]     6.80
##    y0[20]    17.07
##    y0[21]    90.07
##    y0[22]     6.58
##    y0[23]     7.32
##    y0[24]    16.30
##    y0[25]    20.15
##    y0[26]    30.78
##    y0[27]    30.57
##    y0[28]     1.58
##    y0[29]    10.41
##    y0[30]    13.10
##    y0[31]    23.98
##    y0[32]    20.03
##    y0[33]    14.10
##    y0[34]    10.48
##    y0[35]    15.15
##    y0[36]    29.29
##    y0[37]    22.81
##    y0[38]     9.05
##    y0[39]    33.38
##    y0[40]    41.46
##    y0[41]     8.26
##    y0[42]    35.52
##    y0[43]    36.35
##    y0[44]    51.51
##    y0[45]    50.71
##    y0[46]    25.54
##    y0[47]    22.70
##    y0[48]    45.46
##    y0[49]    23.30
##    y0[50]    33.72
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 1.3 seconds.
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 1.3 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 1.4 seconds.
## Chain 4 finished in 1.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.6 seconds.
seeMCMC(mcmc,'m')
##  variable    mean  median   sd  mad      q5    q95 rhat ess_bulk ess_tail
##    lp__   -104.71 -104.47 3.79 3.81 -111.21 -99.12 1.00      699     1158
##    b0[1]    11.40   11.40 6.60 6.67    0.44  21.88 1.00     1499     1037
##    b0[2]     6.48    6.57 4.21 4.34   -0.47  13.22 1.01     2554     1445
##    b0[3]     3.40    3.23 7.75 7.39   -9.48  15.99 1.00     1207     1207
##    b0[4]    10.58   10.49 5.51 5.62    1.63  19.64 1.00     1822     1306
##    b0[5]     8.80    8.80 3.67 3.58    2.78  14.97 1.00     3405     1698
##    b0[6]    14.10   14.10 5.86 5.55    4.30  23.84 1.00     1572     1327
##    b0[7]    20.84   20.87 2.87 2.68   15.90  25.64 1.00     3547     1304
##    b0[8]     5.89    6.02 4.49 4.79   -1.65  12.87 1.00     2184     1242
##    b0[9]    12.59   12.59 4.52 4.49    5.29  20.14 1.00     2125     1367
##    b1[1]     0.41    0.40 1.28 1.25   -1.62   2.52 1.00     1599     1310
##    b1[2]     2.28    2.27 0.44 0.43    1.57   3.02 1.00     2449     1728
##    b1[3]     2.17    2.18 0.69 0.65    1.04   3.34 1.00     1276     1249
##    b1[4]     1.95    1.95 0.41 0.41    1.26   2.62 1.00     1855     1681
##    b1[5]     2.07    2.06 0.40 0.40    1.43   2.72 1.00     2655     1881
##    b1[6]     0.22    0.23 0.46 0.44   -0.55   0.99 1.00     1570     1310
##    b1[7]     3.99    3.99 0.39 0.38    3.37   4.63 1.00     2796     1691
##    b1[8]     1.49    1.48 0.39 0.40    0.89   2.11 1.00     2171     1616
##    b1[9]    -0.67   -0.67 0.37 0.37   -1.26  -0.07 1.00     2019     1596
##    s         5.17    5.10 0.65 0.64    4.19   6.33 1.00     1098     1636
##    m0[1]    16.48   16.40 2.82 2.76   11.82  21.13 1.00     3148     1573
##    m0[2]    42.67   42.68 1.81 1.74   39.68  45.57 1.00     2572     1497
##    m0[3]    36.18   36.16 2.65 2.69   31.78  40.58 1.00     2035     1396
##    m0[4]    12.62   12.52 5.15 4.99    4.16  21.14 1.00     1234     1281
##    m0[5]    10.22   10.22 3.40 3.40    4.85  15.97 1.00     2235     1341
##    m0[6]    27.52   27.48 2.07 2.10   24.18  30.91 1.00     2156     1731
##    m0[7]    13.10   12.98 5.04 5.15    5.02  21.31 1.00     1824     1365
##    m0[8]     3.15    3.20 2.29 2.22   -0.57   6.77 1.00     1953     1173
##    m0[9]    15.37   15.37 3.41 3.25    9.52  21.04 1.00     1627     1404
##    m0[10]   18.44   18.39 3.98 3.79   11.86  25.06 1.00     1700     1420
##    m0[11]    7.08    7.16 4.11 4.24    0.31  13.64 1.00     2553     1446
##    m0[12]   30.03   30.01 3.17 3.22   24.84  35.33 1.00     1788     1449
##    m0[13]   38.19   38.20 1.89 1.84   34.97  41.21 1.00     2855     1622
##    m0[14]    2.43    2.50 2.51 2.46   -1.71   6.45 1.00     1915     1192
##    m0[15]    9.01    9.03 2.88 2.86    4.45  13.90 1.00     2207     1442
##    m0[16]   17.00   17.02 1.87 1.86   13.96  19.96 1.01     2087     1396
##    m0[17]   43.72   43.73 3.38 3.46   38.09  49.15 1.00     2008     1386
##    m0[18]   27.96   27.92 2.08 2.12   24.60  31.35 1.00     2145     1735
##    m0[19]   16.44   16.42 1.90 1.89   13.31  19.48 1.00     1944     1591
##    m0[20]   17.92   17.94 2.28 2.26   14.08  21.62 1.00     2112     1734
##    m0[21]   90.25   90.27 4.85 4.81   82.29  98.45 1.00     2268     1606
##    m0[22]   -0.77   -0.83 3.87 3.78   -7.34   5.53 1.00     1900     1317
##    m0[23]    6.23    6.23 2.06 1.97    2.77   9.60 1.00     2103     1648
##    m0[24]   14.21   14.16 3.02 2.95    9.18  19.24 1.00     3291     1435
##    m0[25]   16.80   16.83 1.79 1.77   13.88  19.72 1.01     2098     1443
##    m0[26]   28.98   29.00 2.31 2.19   25.16  32.74 1.00     3379     1375
##    m0[27]   34.78   34.76 4.02 4.03   28.29  41.63 1.00     1695     1340
##    m0[28]    6.58    6.72 4.34 4.60   -0.71  13.28 1.00     2178     1241
##    m0[29]   12.39   12.35 3.22 3.13    7.05  17.74 1.00     3357     1561
##    m0[30]   18.46   18.45 2.47 2.45   14.40  22.56 1.00     2427     1520
##    m0[31]   26.36   26.36 2.04 2.08   23.03  29.73 1.00     2177     1641
##    m0[32]   19.20   19.21 2.18 2.19   15.54  22.70 1.00     2116     1865
##    m0[33]   12.79   12.85 3.46 3.47    7.17  18.54 1.00     1754     1325
##    m0[34]   17.46   17.48 2.34 2.28   13.55  21.25 1.01     1914     1364
##    m0[35]   16.19   16.23 2.14 2.09   12.61  19.71 1.00     1811     1567
##    m0[36]   32.50   32.53 2.11 2.06   28.96  35.82 1.00     3208     1392
##    m0[37]   23.55   23.46 2.19 2.13   19.94  27.22 1.00     2132     1621
##    m0[38]   12.35   12.43 4.23 4.23    5.19  19.20 1.00     1581     1200
##    m0[39]   31.54   31.57 2.16 2.10   27.95  34.95 1.00     3248     1372
##    m0[40]   42.00   42.06 3.16 3.16   36.67  47.11 1.00     2019     1405
##    m0[41]    6.44    6.41 2.09 2.02    2.92   9.93 1.00     2114     1583
##    m0[42]   28.58   28.52 2.84 2.78   24.15  33.26 1.00     2193     1926
##    m0[43]   32.18   32.15 3.53 3.52   26.62  37.91 1.00     2249     1720
##    m0[44]   45.28   45.38 5.09 5.16   37.02  53.47 1.00     2188     1447
##    m0[45]   43.19   43.20 1.81 1.74   40.19  46.10 1.00     2537     1474
##    m0[46]   26.76   26.73 2.05 2.08   23.42  30.10 1.00     2168     1707
##    m0[47]   16.42   16.42 1.91 1.89   13.28  19.48 1.00     1938     1647
##    m0[48]   44.98   44.99 4.21 4.03   38.44  52.02 1.00     2223     1653
##    m0[49]   14.56   14.56 5.25 5.00    6.09  23.04 1.00     2090     1407
##    m0[50]   37.92   37.92 1.89 1.84   34.70  40.94 1.00     2874     1623
##    y0[1]    16.60   16.56 5.90 5.52    7.18  26.59 1.00     2227     1974
##    y0[2]    42.68   42.60 5.59 5.61   33.52  51.99 1.00     1914     1894
##    y0[3]    36.23   36.24 5.88 5.49   26.42  46.09 1.00     2186     1958
##    y0[4]    12.38   12.38 7.34 7.31    0.24  24.32 1.00     1438     1519
##    y0[5]    10.16   10.16 6.16 6.15    0.06  19.88 1.00     1955     1840
##    y0[6]    27.48   27.67 5.66 5.64   18.07  36.59 1.00     2003     1857
##    y0[7]    13.01   13.02 7.25 7.08    1.26  25.04 1.00     1948     1918
##    y0[8]     3.20    3.12 5.76 5.87   -6.49  12.73 1.00     2035     1768
##    y0[9]    15.16   15.20 6.26 6.26    5.23  25.35 1.00     1888     1886
##    y0[10]   18.39   18.32 6.55 6.59    7.49  29.71 1.00     1753     1776
##    y0[11]    6.99    7.14 6.45 6.46   -3.80  17.46 1.00     1834     1742
##    y0[12]   30.11   30.04 6.09 6.07   20.32  40.03 1.00     2107     1895
##    y0[13]   38.20   38.17 5.59 5.48   29.11  47.37 1.00     2149     1994
##    y0[14]    2.29    2.34 5.86 5.76   -7.47  11.98 1.00     1967     1881
##    y0[15]    8.92    8.81 5.79 5.46   -0.49  18.34 1.00     2062     1658
##    y0[16]   16.96   16.95 5.54 5.40    7.75  25.97 1.00     1879     1547
##    y0[17]   43.86   43.76 6.27 6.02   33.59  54.19 1.00     1862     2013
##    y0[18]   27.98   27.93 5.55 5.40   19.10  37.11 1.00     2018     1959
##    y0[19]   16.42   16.24 5.47 5.39    7.55  25.62 1.00     2114     2012
##    y0[20]   17.94   17.87 5.48 5.29    8.94  27.29 1.00     1964     1732
##    y0[21]   90.23   90.05 7.03 6.95   78.71 101.87 1.00     2317     2153
##    y0[22]   -0.71   -0.65 6.46 6.37  -11.24  10.11 1.00     2042     1878
##    y0[23]    6.23    6.20 5.62 5.65   -3.18  15.36 1.00     2077     1440
##    y0[24]   14.11   14.04 6.19 6.04    3.79  24.11 1.00     2065     1996
##    y0[25]   16.77   16.72 5.47 5.46    7.67  25.57 1.00     1780     1825
##    y0[26]   28.88   28.87 5.76 5.68   19.37  38.40 1.00     2047     1962
##    y0[27]   34.66   34.63 6.48 6.35   23.86  45.49 1.00     1884     1853
##    y0[28]    6.71    6.75 6.64 6.41   -4.37  17.69 1.00     2125     1972
##    y0[29]   12.30   12.25 6.20 6.22    2.23  22.14 1.00     2118     1949
##    y0[30]   18.51   18.54 5.80 5.79    9.20  28.21 1.00     2089     1881
##    y0[31]   26.24   26.13 5.86 6.04   16.58  35.83 1.00     1997     1805
##    y0[32]   19.24   19.29 5.66 5.32    9.96  28.60 1.00     2002     1970
##    y0[33]   12.80   12.84 6.29 6.32    2.77  23.03 1.00     2046     1986
##    y0[34]   17.42   17.36 5.58 5.47    8.34  26.70 1.00     2097     1913
##    y0[35]   16.23   16.16 5.80 5.80    6.76  25.73 1.00     1963     1933
##    y0[36]   32.47   32.43 5.36 5.29   23.85  41.55 1.00     2091     2042
##    y0[37]   23.53   23.54 5.67 5.46   13.89  32.71 1.00     1880     1855
##    y0[38]   12.25   12.35 6.65 6.69    1.14  23.07 1.00     1709     1741
##    y0[39]   31.54   31.60 5.63 5.72   22.22  40.51 1.00     2305     2012
##    y0[40]   42.09   42.33 6.15 6.15   31.94  52.08 1.00     2033     1924
##    y0[41]    6.12    6.19 5.60 5.44   -3.16  15.26 1.00     2161     2015
##    y0[42]   28.60   28.47 5.98 5.94   19.15  38.09 1.01     1883     1881
##    y0[43]   32.09   32.02 6.13 6.10   22.13  42.24 1.00     2303     2001
##    y0[44]   45.27   45.39 7.38 7.23   32.60  57.25 1.00     1996     1704
##    y0[45]   43.15   43.13 5.41 5.37   34.54  51.99 1.00     1973     1973
##    y0[46]   26.64   26.67 5.56 5.47   17.50  35.82 1.00     2013     1893
##    y0[47]   16.34   16.38 5.54 5.49    7.40  25.44 1.00     1887     1819
##    y0[48]   44.94   44.83 6.61 6.42   33.83  56.10 1.00     2295     1912
##    y0[49]   14.47   14.25 7.72 7.49    2.06  27.44 1.00     1981     1416
##    y0[50]   37.80   37.74 5.47 5.23   28.55  46.76 1.00     2121     2151

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

estimate as class have relation each other  
all b0l,b1l have a distribution  
class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)

ex12-2.stan

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  real b00;
  real<lower=0,upper=100> sb0;
  vector[K] b0;
  real b10;
  real<lower=0,upper=100> sb1;
  vector[K] b1;
  real<lower=0,upper=100> s;
}
model{
  b0~normal(b00,sb0);
  b1~normal(b10,sb1);
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-2.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -266.657 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -39.1515     0.0339647        538871      0.2692      0.2692      150    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      122      -33.8812   0.000129671       547.855           1           1      178    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
try(print(mle))
##  variable estimate
##    lp__     -33.88
##    b00        1.17
##    sb0        0.00
##    b0[1]      1.17
##    b0[2]      1.17
##    b0[3]      1.17
##    b0[4]      1.17
##    b0[5]      1.17
##    b0[6]      1.17
##    b0[7]      1.17
##    b0[8]      1.17
##    b0[9]      1.17
##    b10        2.50
##    sb1        0.45
##    b1[1]      4.14
##    b1[2]      2.94
##    b1[3]      2.60
##    b1[4]      2.52
##    b1[5]      2.09
##    b1[6]      1.57
##    b1[7]      5.53
##    b1[8]      1.70
##    b1[9]      0.29
##    s         11.44
##    m0[1]      8.91
##    m0[2]     31.38
##    m0[3]     34.27
##    m0[4]     12.22
##    m0[5]      2.18
##    m0[6]     28.36
##    m0[7]      4.43
##    m0[8]      5.21
##    m0[9]     10.24
##    m0[10]    32.07
##    m0[11]     1.95
##    m0[12]    33.08
##    m0[13]    25.19
##    m0[14]     5.51
##    m0[15]     2.70
##    m0[16]    21.84
##    m0[17]    44.03
##    m0[18]    28.93
##    m0[19]    17.83
##    m0[20]    14.92
##    m0[21]    97.23
##    m0[22]     6.89
##    m0[23]     3.89
##    m0[24]     6.62
##    m0[25]    20.40
##    m0[26]    12.44
##    m0[27]    38.76
##    m0[28]     1.96
##    m0[29]     4.79
##    m0[30]    16.66
##    m0[31]    26.86
##    m0[32]    16.38
##    m0[33]    15.44
##    m0[34]    25.15
##    m0[35]    16.09
##    m0[36]    17.31
##    m0[37]    21.35
##    m0[38]    10.88
##    m0[39]    15.98
##    m0[40]    41.80
##    m0[41]     3.80
##    m0[42]    27.09
##    m0[43]    31.21
##    m0[44]    37.92
##    m0[45]    32.11
##    m0[46]    27.38
##    m0[47]    17.75
##    m0[48]    50.92
##    m0[49]    33.55
##    m0[50]    24.81
##    y0[1]     16.40
##    y0[2]     22.61
##    y0[3]     43.99
##    y0[4]     -4.37
##    y0[5]      5.40
##    y0[6]     11.56
##    y0[7]      9.39
##    y0[8]     -7.71
##    y0[9]     11.41
##    y0[10]    51.98
##    y0[11]     3.20
##    y0[12]    11.05
##    y0[13]    51.59
##    y0[14]    13.30
##    y0[15]   -11.74
##    y0[16]     7.34
##    y0[17]    21.90
##    y0[18]    45.90
##    y0[19]    -6.45
##    y0[20]     4.72
##    y0[21]   108.35
##    y0[22]    15.75
##    y0[23]   -20.95
##    y0[24]    -0.88
##    y0[25]    32.32
##    y0[26]     1.66
##    y0[27]    54.45
##    y0[28]     8.24
##    y0[29]     0.68
##    y0[30]    32.66
##    y0[31]    47.53
##    y0[32]    14.66
##    y0[33]    20.98
##    y0[34]    24.57
##    y0[35]    19.70
##    y0[36]    21.27
##    y0[37]    10.48
##    y0[38]     3.96
##    y0[39]    28.88
##    y0[40]    31.95
##    y0[41]    -5.74
##    y0[42]    37.84
##    y0[43]    27.02
##    y0[44]    53.68
##    y0[45]    28.46
##    y0[46]    30.82
##    y0[47]    -2.94
##    y0[48]    49.67
##    y0[49]    28.61
##    y0[50]    39.75
mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 1 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 2 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 3 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 3 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 4 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 4 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 4 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 3 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 3 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 4 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 1 finished in 1.2 seconds.
## Chain 3 finished in 1.2 seconds.
## Chain 4 finished in 1.1 seconds.
## Chain 2 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 2 finished in 1.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.2 seconds.
## Total execution time: 1.4 seconds.
seeMCMC(mcmc,'m')
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -128.16 -127.75 4.45 4.36 -136.09 -121.56 1.00     1765     2672
##    b00      10.85   10.93 2.65 2.45    6.45   14.92 1.00     4954     5420
##    sb0       5.53    5.12 2.54 2.11    2.13   10.16 1.00     2183     1828
##    b0[1]    10.22   10.29 3.97 3.77    3.59   16.58 1.00     6858     5653
##    b0[2]     8.52    8.60 3.29 3.25    3.05   13.78 1.00     5341     5107
##    b0[3]     8.52    8.76 4.59 4.37    0.70   15.47 1.00     5075     4803
##    b0[4]    10.88   10.96 3.82 3.60    4.49   17.11 1.00     6748     5749
##    b0[5]     9.75    9.76 2.99 2.91    4.83   14.58 1.00     6680     6133
##    b0[6]    11.69   11.71 4.12 3.87    5.02   18.54 1.00     5907     5175
##    b0[7]    18.93   18.98 3.05 3.03   13.80   23.84 1.00     3737     3172
##    b0[8]     8.09    8.25 3.65 3.53    1.93   13.72 1.00     4230     5431
##    b0[9]    11.11   11.18 3.46 3.34    5.29   16.68 1.00     6632     5513
##    b10       1.54    1.54 0.61 0.54    0.56    2.51 1.00     6480     4133
##    sb1       1.66    1.54 0.58 0.46    0.99    2.71 1.00     6736     4916
##    b1[1]     0.74    0.73 0.82 0.81   -0.61    2.07 1.00     6883     5807
##    b1[2]     2.07    2.07 0.36 0.35    1.49    2.68 1.00     5888     5471
##    b1[3]     1.74    1.73 0.45 0.44    1.04    2.48 1.00     5833     4721
##    b1[4]     1.93    1.93 0.31 0.30    1.42    2.44 1.00     7114     6069
##    b1[5]     2.00    2.00 0.35 0.35    1.41    2.58 1.00     7688     5525
##    b1[6]     0.41    0.41 0.34 0.33   -0.14    0.97 1.00     6121     5291
##    b1[7]     4.13    4.13 0.41 0.40    3.48    4.82 1.00     4370     3119
##    b1[8]     1.33    1.33 0.32 0.32    0.82    1.87 1.00     5195     5527
##    b1[9]    -0.53   -0.54 0.30 0.30   -1.02   -0.03 1.00     7064     5650
##    s         5.06    5.00 0.62 0.61    4.15    6.17 1.00     5593     5711
##    m0[1]    17.16   17.16 2.40 2.35   13.24   21.03 1.00     7478     6456
##    m0[2]    41.54   41.57 1.85 1.83   38.40   44.52 1.00     6204     5545
##    m0[3]    36.17   36.16 2.51 2.44   32.10   40.32 1.00     9774     5518
##    m0[4]    15.93   16.06 3.22 3.12   10.43   21.00 1.00     5051     5155
##    m0[5]     9.22    9.27 2.62 2.54    4.83   13.46 1.00     6268     6157
##    m0[6]    27.71   27.68 1.92 1.88   24.57   30.90 1.00     7107     5251
##    m0[7]    13.37   13.44 3.50 3.31    7.55   19.07 1.00     6807     5928
##    m0[8]     3.61    3.61 2.25 2.20   -0.07    7.34 1.00     8597     6083
##    m0[9]    14.09   14.06 2.48 2.34   10.05   18.24 1.00     6139     5955
##    m0[10]   19.86   19.86 3.38 3.33   14.30   25.40 1.00     7080     5864
##    m0[11]    9.07    9.14 3.22 3.18    3.72   14.21 1.00     5341     5111
##    m0[12]   29.91   29.88 3.11 3.09   24.86   35.08 1.00     7795     5826
##    m0[13]   36.90   36.94 1.95 1.93   33.58   40.05 1.00     4972     4535
##    m0[14]    3.04    3.03 2.45 2.43   -0.97    7.11 1.00     8506     6209
##    m0[15]    8.26    8.30 2.27 2.22    4.42   11.91 1.00     6388     6281
##    m0[16]   17.16   17.16 1.86 1.86   14.11   20.19 1.00     8137     5637
##    m0[17]   43.63   43.64 3.20 3.19   38.36   48.90 1.00     9451     5870
##    m0[18]   28.11   28.09 1.95 1.91   24.92   31.34 1.00     7157     5139
##    m0[19]   16.10   16.09 1.73 1.70   13.23   18.89 1.00     7500     5043
##    m0[20]   18.84   18.85 2.03 2.00   15.45   22.17 1.00     6902     5863
##    m0[21]   90.81   90.82 5.09 4.95   82.39   99.00 1.00     5772     5411
##    m0[22]    0.50    0.44 3.56 3.53   -5.27    6.32 1.00     8047     6258
##    m0[23]    6.06    6.06 1.84 1.80    3.02    9.14 1.00     8356     5358
##    m0[24]   14.97   14.99 2.52 2.46   10.85   19.04 1.00     7126     6251
##    m0[25]   16.77   16.78 1.77 1.77   13.88   19.64 1.00     8012     5554
##    m0[26]   27.37   27.44 2.44 2.39   23.21   31.26 1.00     3957     3303
##    m0[27]   33.72   33.70 3.73 3.64   27.56   39.89 1.00     7685     5999
##    m0[28]    8.71    8.86 3.52 3.42    2.73   14.14 1.00     4262     5560
##    m0[29]   13.22   13.23 2.65 2.60    8.87   17.47 1.00     6931     6220
##    m0[30]   19.45   19.46 2.02 1.99   16.09   22.71 1.00     5750     6319
##    m0[31]   26.65   26.62 1.88 1.83   23.57   29.74 1.00     6951     5489
##    m0[32]   19.98   19.99 1.99 1.96   16.73   23.26 1.00     7441     5842
##    m0[33]   12.76   12.76 2.72 2.67    8.20   17.19 1.00     8731     6005
##    m0[34]   18.03   18.03 2.23 2.22   14.39   21.65 1.00     7863     6077
##    m0[35]   15.63   15.63 1.80 1.76   12.73   18.58 1.00     7042     5930
##    m0[36]   31.01   31.07 2.22 2.18   27.23   34.53 1.00     4200     3636
##    m0[37]   23.86   23.87 2.12 2.10   20.44   27.42 1.00     8761     5481
##    m0[38]   11.94   11.98 2.89 2.77    7.04   16.59 1.00     7887     6598
##    m0[39]   30.02   30.08 2.27 2.24   26.15   33.63 1.00     4123     3702
##    m0[40]   41.93   41.93 3.01 3.02   36.97   46.90 1.00     9575     5830
##    m0[41]    6.23    6.24 1.85 1.80    3.18    9.30 1.00     8288     5600
##    m0[42]   28.35   28.35 2.71 2.69   23.86   32.82 1.00     8533     6366
##    m0[43]   31.56   31.57 3.29 3.29   26.15   36.96 1.00     8047     6589
##    m0[44]   44.95   44.93 4.94 4.80   36.86   53.08 1.00     9177     6502
##    m0[45]   42.08   42.12 1.85 1.82   38.95   45.07 1.00     6381     5525
##    m0[46]   27.01   26.98 1.89 1.85   23.90   30.14 1.00     7009     5374
##    m0[47]   16.07   16.07 1.73 1.70   13.21   18.87 1.00     7472     5019
##    m0[48]   43.62   43.63 3.81 3.77   37.39   49.94 1.00     6826     5018
##    m0[49]   15.98   15.98 4.44 4.41    8.75   23.35 1.00     8796     6041
##    m0[50]   36.62   36.65 1.96 1.93   33.28   39.77 1.00     4921     4562
##    y0[1]    17.17   17.20 5.69 5.64    7.71   26.36 1.00     7950     7366
##    y0[2]    41.54   41.58 5.44 5.33   32.52   50.25 1.00     8507     7925
##    y0[3]    36.19   36.17 5.70 5.58   26.86   45.53 1.00     8484     7649
##    y0[4]    16.00   16.06 6.08 5.91    5.79   25.84 1.00     6841     7117
##    y0[5]     9.14    9.13 5.80 5.75   -0.20   18.75 1.00     7931     7594
##    y0[6]    27.73   27.72 5.40 5.22   18.78   36.62 1.00     8138     7539
##    y0[7]    13.37   13.40 6.13 6.01    3.27   23.41 1.00     7638     7251
##    y0[8]     3.56    3.54 5.62 5.47   -5.51   12.83 1.00     8030     7906
##    y0[9]    14.21   14.16 5.63 5.59    4.84   23.49 1.00     7756     7567
##    y0[10]   19.89   19.88 6.11 6.03    9.80   29.76 1.00     7775     6792
##    y0[11]    9.11    9.18 6.03 6.00   -0.76   19.03 1.00     6688     7255
##    y0[12]   29.86   29.92 5.95 5.82   20.10   39.66 1.00     7890     7592
##    y0[13]   36.97   37.06 5.38 5.31   28.03   45.78 1.00     7755     7337
##    y0[14]    3.06    3.03 5.66 5.51   -6.12   12.39 1.00     8069     7334
##    y0[15]    8.32    8.29 5.51 5.46   -0.70   17.44 1.00     7733     7628
##    y0[16]   17.15   17.08 5.43 5.43    8.37   26.01 1.00     7932     7557
##    y0[17]   43.62   43.57 6.06 5.92   33.62   53.56 1.00     8761     7604
##    y0[18]   28.06   28.03 5.54 5.36   19.29   37.33 1.00     8136     7810
##    y0[19]   16.06   16.08 5.29 5.20    7.38   24.78 1.00     8001     7505
##    y0[20]   18.81   18.82 5.46 5.38    9.87   27.78 1.00     7936     7910
##    y0[21]   90.76   90.76 7.20 6.97   79.01  102.75 1.00     7013     7182
##    y0[22]    0.53    0.56 6.19 6.07   -9.68   10.73 1.00     8004     7789
##    y0[23]    6.03    5.99 5.37 5.27   -2.78   14.97 1.00     7943     7643
##    y0[24]   14.99   14.92 5.64 5.66    5.78   24.21 1.00     8510     7515
##    y0[25]   16.80   16.74 5.36 5.21    8.02   25.68 1.00     8028     7809
##    y0[26]   27.44   27.40 5.70 5.56   18.01   36.96 1.00     7021     6967
##    y0[27]   33.66   33.69 6.34 6.24   23.39   44.21 1.00     8094     7829
##    y0[28]    8.59    8.60 6.24 6.13   -1.66   18.74 1.00     7124     6994
##    y0[29]   13.17   13.19 5.67 5.62    3.83   22.31 1.00     7535     7267
##    y0[30]   19.43   19.39 5.45 5.35   10.70   28.37 1.00     7289     7584
##    y0[31]   26.67   26.68 5.41 5.31   17.83   35.57 1.00     8384     7017
##    y0[32]   19.97   19.96 5.48 5.36   11.02   28.97 1.00     8102     7887
##    y0[33]   12.74   12.76 5.81 5.71    3.34   22.35 1.00     8199     7693
##    y0[34]   17.98   18.00 5.61 5.54    8.81   27.38 1.00     8084     7740
##    y0[35]   15.68   15.69 5.42 5.32    6.76   24.49 1.00     8160     7920
##    y0[36]   31.04   31.10 5.58 5.48   21.91   40.22 1.00     7240     7178
##    y0[37]   23.86   23.84 5.42 5.36   14.87   32.70 1.00     7803     7716
##    y0[38]   11.97   12.01 5.92 5.82    2.30   21.64 1.00     8236     8053
##    y0[39]   30.03   30.07 5.55 5.41   20.95   39.16 1.00     7440     7646
##    y0[40]   42.02   42.10 5.92 5.75   32.26   51.83 1.00     8188     7657
##    y0[41]    6.30    6.29 5.49 5.40   -2.80   15.47 1.00     7765     7888
##    y0[42]   28.35   28.36 5.74 5.68   19.09   37.60 1.00     8085     8013
##    y0[43]   31.64   31.70 5.99 5.93   21.51   41.34 1.00     7802     7412
##    y0[44]   44.96   44.86 7.18 7.15   33.18   56.70 1.00     8643     6990
##    y0[45]   41.99   41.95 5.45 5.32   33.06   51.00 1.00     7661     7369
##    y0[46]   27.05   27.03 5.45 5.38   18.13   35.84 1.00     7599     7636
##    y0[47]   16.03   16.07 5.30 5.12    7.27   24.59 1.00     7817     7339
##    y0[48]   43.54   43.58 6.37 6.39   32.91   53.83 1.00     8179     7815
##    y0[49]   16.11   16.12 6.77 6.68    4.87   27.17 1.00     8543     7976
##    y0[50]   36.55   36.65 5.48 5.33   27.44   45.47 1.00     6647     7441

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

ex13

generalized linear mixed model

(X,y)i=1-n
b[b0,b1,...]
linear model    y~N(Xb,s)  
generalized linear model    y~dist.(m=link(Xb),s)  

fixed effect    b0, b1,...  
individual random effect   b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...  
class c=1-k  
class effect    b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...  

for y=dist.(m,s)
random intercept model    m=b0+r0i+b1*x, m=b0+r0c+b1*x  
random coefficient model  m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x  
mixed model   m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x  

note  
@ yi=b0+b1*xi+r0i is not useful, ri is included in s  
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
  but variance is larger than mean (over dispersion),
  random intercept model is useful
n=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)

ex13-0.stan

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
}
model{
  y~poisson_log(b0+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-0.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -3304.39 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12       658.086   0.000343231     0.0028631           1           1       17    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     658.09
##    b0         2.17
##    b1         0.11
##    m0[1]      2.22
##    m0[2]      2.69
##    m0[3]      2.22
##    m0[4]      2.25
##    m0[5]      3.00
##    m0[6]      3.30
##    m0[7]      2.88
##    m0[8]      2.68
##    m0[9]      3.14
##    m0[10]     3.00
##    m0[11]     2.88
##    m0[12]     3.03
##    m0[13]     3.24
##    m0[14]     2.39
##    m0[15]     2.37
##    m0[16]     3.20
##    m0[17]     2.77
##    m0[18]     3.30
##    m0[19]     2.62
##    m0[20]     2.51
##    y0[1]      9.00
##    y0[2]     13.00
##    y0[3]     14.00
##    y0[4]      5.00
##    y0[5]     18.00
##    y0[6]     30.00
##    y0[7]     20.00
##    y0[8]     16.00
##    y0[9]     26.00
##    y0[10]    21.00
##    y0[11]    21.00
##    y0[12]    24.00
##    y0[13]    30.00
##    y0[14]    16.00
##    y0[15]     8.00
##    y0[16]    27.00
##    y0[17]    15.00
##    y0[18]    21.00
##    y0[19]    19.00
##    y0[20]     5.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   657.13 657.44 1.00 0.67 655.08 658.04 1.01      652      574
##    b0       2.17   2.17 0.12 0.12   1.97   2.37 1.01      494      394
##    b1       0.11   0.11 0.02 0.02   0.09   0.14 1.01      604      619
##    m0[1]    2.22   2.22 0.12 0.12   2.03   2.41 1.01      494      388
##    m0[2]    2.69   2.69 0.06 0.06   2.59   2.79 1.01      638      611
##    m0[3]    2.22   2.23 0.12 0.12   2.03   2.41 1.01      494      387
##    m0[4]    2.25   2.25 0.11 0.11   2.06   2.43 1.01      493      388
##    m0[5]    2.99   3.00 0.06 0.05   2.90   3.09 1.00     1886     1332
##    m0[6]    3.30   3.30 0.08 0.08   3.16   3.43 1.00     1371     1383
##    m0[7]    2.88   2.88 0.05 0.05   2.79   2.97 1.00     1283     1363
##    m0[8]    2.67   2.67 0.06 0.06   2.57   2.78 1.01      616      512
##    m0[9]    3.13   3.14 0.07 0.06   3.02   3.24 1.00     1828     1466
##    m0[10]   3.00   3.00 0.06 0.05   2.90   3.09 1.00     1897     1332
##    m0[11]   2.88   2.88 0.05 0.05   2.79   2.97 1.00     1293     1363
##    m0[12]   3.03   3.03 0.06 0.06   2.93   3.12 1.00     1972     1633
##    m0[13]   3.24   3.24 0.08 0.07   3.11   3.36 1.00     1516     1424
##    m0[14]   2.39   2.39 0.09 0.09   2.24   2.55 1.01      495      365
##    m0[15]   2.37   2.37 0.10 0.09   2.21   2.53 1.01      494      370
##    m0[16]   3.19   3.20 0.07 0.07   3.07   3.30 1.00     1657     1445
##    m0[17]   2.77   2.77 0.06 0.06   2.67   2.86 1.01      786      764
##    m0[18]   3.30   3.30 0.08 0.08   3.16   3.43 1.00     1365     1383
##    m0[19]   2.62   2.62 0.07 0.07   2.50   2.73 1.01      561      483
##    m0[20]   2.51   2.51 0.08 0.08   2.38   2.64 1.01      511      393
##    y0[1]    9.19   9.00 3.15 2.97   4.00  15.00 1.00     1746     1725
##    y0[2]   14.73  15.00 3.94 4.45   9.00  21.00 1.00     1923     2035
##    y0[3]    9.26   9.00 3.23 2.97   4.00  15.00 1.00     1612     1931
##    y0[4]    9.54   9.00 3.30 2.97   5.00  16.00 1.00     1426     1499
##    y0[5]   19.92  20.00 4.69 4.45  13.00  28.00 1.00     1929     1807
##    y0[6]   27.22  27.00 5.64 5.93  18.00  37.00 1.00     2115     2037
##    y0[7]   17.80  18.00 4.38 4.45  11.00  25.00 1.00     1948     2025
##    y0[8]   14.57  14.00 4.00 4.45   8.00  22.00 1.00     1825     2074
##    y0[9]   22.90  23.00 4.87 4.45  15.00  31.00 1.00     1764     1969
##    y0[10]  20.00  20.00 4.54 4.45  13.00  28.00 1.00     1969     1879
##    y0[11]  17.85  18.00 4.34 4.45  11.00  25.00 1.00     1888     1972
##    y0[12]  20.76  20.00 4.75 4.45  13.00  29.00 1.00     2076     1848
##    y0[13]  25.79  26.00 5.35 5.93  17.00  35.00 1.00     1922     1837
##    y0[14]  11.03  11.00 3.45 2.97   6.00  17.00 1.00     1682     1927
##    y0[15]  10.62  10.00 3.31 2.97   6.00  16.00 1.00     1604     1925
##    y0[16]  24.42  24.00 5.38 5.93  16.00  33.00 1.00     1916     1892
##    y0[17]  15.96  16.00 4.18 4.45   9.95  23.00 1.00     1915     2002
##    y0[18]  27.17  27.00 5.91 5.93  18.00  37.00 1.00     1537     1929
##    y0[19]  13.73  14.00 3.76 4.45   8.00  20.00 1.00     1556     1721
##    y0[20]  12.49  12.00 3.65 2.97   7.00  19.00 1.00     1547     1854

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

ex13-1.stan

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
  real<lower=0> sr0;
  vector[N] r0;
}
model{
  r0~normal(0,sr0);
  for(i in 1:N)
    y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+r0[i]+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-1.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -77904 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99       13247.1     0.0196361       261.206      0.4196      0.4196      141    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      148       13353.4    0.00102166          6298      0.5652      0.5652      221    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__   13353.40
##    b0         1.73
##    b1         0.17
##    sr0        0.00
##    r0[1]      0.00
##    r0[2]      0.00
##    r0[3]      0.00
##    r0[4]      0.00
##    r0[5]      0.00
##    r0[6]      0.00
##    r0[7]      0.00
##    r0[8]      0.00
##    r0[9]      0.00
##    r0[10]     0.00
##    r0[11]     0.00
##    r0[12]     0.00
##    r0[13]     0.00
##    r0[14]     0.00
##    r0[15]     0.00
##    r0[16]     0.00
##    r0[17]     0.00
##    r0[18]     0.00
##    r0[19]     0.00
##    r0[20]     0.00
##    m0[1]      1.81
##    m0[2]      2.52
##    m0[3]      1.81
##    m0[4]      1.84
##    m0[5]      2.98
##    m0[6]      3.45
##    m0[7]      2.81
##    m0[8]      2.50
##    m0[9]      3.20
##    m0[10]     2.99
##    m0[11]     2.81
##    m0[12]     3.04
##    m0[13]     3.36
##    m0[14]     2.07
##    m0[15]     2.03
##    m0[16]     3.29
##    m0[17]     2.64
##    m0[18]     3.45
##    m0[19]     2.41
##    m0[20]     2.24
##    y0[1]      9.00
##    y0[2]     15.00
##    y0[3]      4.00
##    y0[4]      6.00
##    y0[5]     19.00
##    y0[6]     26.00
##    y0[7]     11.00
##    y0[8]     12.00
##    y0[9]     24.00
##    y0[10]    20.00
##    y0[11]     8.00
##    y0[12]    20.00
##    y0[13]    30.00
##    y0[14]     5.00
##    y0[15]     7.00
##    y0[16]    32.00
##    y0[17]     8.00
##    y0[18]    39.00
##    y0[19]    13.00
##    y0[20]    15.00
mcmc=goMCMC(mdl,data,wrm=500)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.9 seconds.
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 finished in 0.9 seconds.
## Chain 3 finished in 0.9 seconds.
## Chain 4 finished in 0.9 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.9 seconds.
## Total execution time: 1.0 seconds.
seeMCMC(mcmc,'[m,r]')
##  variable     mean   median    sd   mad       q5      q95 rhat ess_bulk
##    lp__   13242.54 13242.20 14.49 15.57 13220.30 13268.10 1.05       53
##    b0         2.17     2.17  0.03  0.03     2.12     2.22 1.01      652
##    b1         0.11     0.11  0.00  0.00     0.11     0.12 1.01      671
##    sr0        0.01     0.01  0.01  0.01     0.00     0.02 1.05       55
##    r0[1]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     2669
##    r0[2]      0.00     0.00  0.01  0.01    -0.02     0.02 1.01     3337
##    r0[3]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     2843
##    r0[4]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3511
##    r0[5]      0.00     0.00  0.01  0.01    -0.02     0.02 1.01     3200
##    r0[6]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     2753
##    r0[7]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3004
##    r0[8]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3743
##    r0[9]      0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3754
##    r0[10]     0.00     0.00  0.01  0.01    -0.02     0.02 1.01     4534
##    r0[11]     0.00     0.00  0.01  0.01    -0.02     0.02 1.01     2822
##    r0[12]     0.00     0.00  0.01  0.01    -0.02     0.02 1.01     2662
##    r0[13]     0.00     0.00  0.01  0.01    -0.02     0.02 1.01     2919
##    r0[14]     0.00     0.00  0.01  0.01    -0.02     0.02 1.01     3858
##    r0[15]     0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3170
##    r0[16]     0.00     0.00  0.01  0.01    -0.02     0.02 1.01     3022
##    r0[17]     0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3075
##    r0[18]     0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3845
##    r0[19]     0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3313
##    r0[20]     0.00     0.00  0.01  0.01    -0.02     0.02 1.02     3432
##    m0[1]      2.22     2.22  0.03  0.03     2.18     2.27 1.01      815
##    m0[2]      2.69     2.69  0.02  0.02     2.66     2.72 1.00     1350
##    m0[3]      2.23     2.23  0.03  0.03     2.18     2.27 1.00      810
##    m0[4]      2.25     2.25  0.03  0.03     2.20     2.29 1.01      797
##    m0[5]      3.00     3.00  0.02  0.01     2.97     3.02 1.00     2272
##    m0[6]      3.30     3.30  0.02  0.02     3.26     3.34 1.00     1098
##    m0[7]      2.88     2.88  0.02  0.01     2.86     2.91 1.00     2001
##    m0[8]      2.68     2.68  0.02  0.02     2.65     2.71 1.00     1267
##    m0[9]      3.14     3.14  0.02  0.02     3.11     3.17 1.00     2138
##    m0[10]     3.00     3.00  0.02  0.02     2.97     3.03 1.00     2745
##    m0[11]     2.88     2.88  0.02  0.02     2.86     2.91 1.00     2114
##    m0[12]     3.03     3.03  0.02  0.02     3.01     3.06 1.00     2245
##    m0[13]     3.24     3.24  0.02  0.02     3.21     3.28 1.00     1705
##    m0[14]     2.39     2.39  0.02  0.03     2.35     2.43 1.01      810
##    m0[15]     2.37     2.37  0.02  0.02     2.33     2.41 1.01      847
##    m0[16]     3.20     3.20  0.02  0.02     3.16     3.23 1.00     1644
##    m0[17]     2.77     2.77  0.02  0.02     2.74     2.80 1.00     1497
##    m0[18]     3.30     3.30  0.02  0.02     3.27     3.34 1.00     1634
##    m0[19]     2.62     2.62  0.02  0.02     2.59     2.65 1.00     1086
##    m0[20]     2.51     2.51  0.02  0.02     2.48     2.54 1.00      906
##    y0[1]      9.28     9.00  2.99  2.97     5.00    14.00 1.00     2112
##    y0[2]     14.73    15.00  3.92  4.45     9.00    22.00 1.00     2117
##    y0[3]      9.23     9.00  3.09  2.97     5.00    15.00 1.00     1902
##    y0[4]      9.49     9.00  3.01  2.97     5.00    15.00 1.00     1986
##    y0[5]     20.06    20.00  4.52  4.45    13.00    28.00 1.00     2115
##    y0[6]     27.27    27.00  5.23  5.93    19.00    36.00 1.00     1904
##    y0[7]     17.82    18.00  4.17  4.45    11.00    25.00 1.00     1979
##    y0[8]     14.57    15.00  3.78  4.45     9.00    21.00 1.00     2096
##    y0[9]     22.90    23.00  4.89  4.45    15.00    31.00 1.00     1798
##    y0[10]    20.20    20.00  4.56  4.45    13.00    28.00 1.00     1876
##    y0[11]    17.87    18.00  4.30  4.45    11.00    25.00 1.00     1946
##    y0[12]    20.88    21.00  4.65  4.45    13.00    29.00 1.00     1852
##    y0[13]    25.59    25.00  5.25  4.45    17.00    34.00 1.00     1675
##    y0[14]    10.98    11.00  3.29  2.97     6.00    17.00 1.00     1874
##    y0[15]    10.73    11.00  3.30  2.97     6.00    16.00 1.00     2061
##    y0[16]    24.45    24.00  5.12  4.45    16.00    33.00 1.00     2019
##    y0[17]    15.77    16.00  3.96  4.45     9.00    23.00 1.00     2034
##    y0[18]    27.10    27.00  5.30  5.93    19.00    36.00 1.00     1937
##    y0[19]    13.75    13.00  3.78  4.45     8.00    20.00 1.00     1851
##    y0[20]    12.36    12.00  3.51  2.97     7.00    18.00 1.00     1943
##  ess_tail
##        87
##      1295
##      1172
##        72
##       901
##       741
##       808
##       793
##       788
##       599
##       839
##       757
##       827
##       829
##       674
##       819
##       763
##       647
##       646
##       719
##       649
##       915
##       762
##       736
##      1418
##      1410
##      1269
##      1324
##      1108
##      1333
##      1304
##      1249
##      1199
##      1345
##      1304
##      1299
##      1336
##      1236
##      1170
##      1305
##       975
##      1364
##      1320
##      1329
##      1918
##      2011
##      1853
##      2089
##      1829
##      1755
##      1758
##      1924
##      1649
##      1948
##      1917
##      1983
##      1854
##      1910
##      1972
##      2051
##      1877
##      1711
##      2001
##      1933

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')